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Abstract 

We revisit the Lagrange and Delaunay systems of equations of celestial mechanics, 
and point out a previously neglected aspect of these equations: in both cases the orbit 
resides on a certain 9(N-l)-dimensional submanifold of the 12(N-l)-dimensional space 
spanned by the orbital elements and their time derivatives. We demonstrate that there 
exists a vast freedom in choosing this submanifold. This freedom of choice (=freedom 
of gauge fixing) reveals a symmetry hiding behind Lagrange's and Delaunay's sys- 
tems, which is, mathematically, analogous to the gauge invariance in electrodynamics. 
Just like a convenient choice of gauge simplifies calculations in electrodynamics, so 
the freedom of choice of the submanifold may, potentially, be used to create simpler 
schemes of orbit integration. On the other hand, the presence of this feature may be 
a previously unrecognised source of numerical instability. 
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1 Introduction 



1.1 The Variation-of-Parameters Method 

The method of variation of parameters (VOP) was invented in the middle of XVIII-th 
century by Euler (1748, 1753) as a tool for treating highly-nonlinear problems emerging 
in the context of celestial mechanics. Though Euler himself eventually lost interest in this 
approach, his line of research was intensively continued by Lagrange who employed it for 
deriving his celebrated system of equations describing evolution of the orbital parameters 
(which in the astronomical parlance are called "orbital elements"). 

In the modern textbooks, the VOP method is normally introduced as one of the means 
by which one can solve an inhomogeneous linear differential equation: one first finds all 
solutions of the appropriate linear homogeneous equation, and then instills time dependence 
into the coefficients in the linear combination of these solutions. Here follows the easiest 
example: 

x + p(t) x + q(t) x = F(t) . (1) 

To solve this inhomogeneous equation, one starts out with the homogeneous one: 

x + p(t) x + q(t) x = . (2) 

A linear combination of its two fundamental solutions will read: S Xi{t) + C X2{t) . The 
recipe has it that at this point one should look for a solution to (jTJ) in ansatz 

x = S(t) Xl (t) + C(t) x 2 (t) . (3) 

Since the functions xi^t) are known, what one has to do is just to find S{t) and C(t). 
Equation (TjQ) will, by itself, be insufficient for determining two independent functions. The 
excessive freedom can be removed through a by-hand imposure of an extra equality often 
chosen as 

Sxj. + C x 2 = . (4) 
It greatly simplifies the expressions for x and x: 

x = S ii + C ±2 , x = S x\ + C ±2 + C x\ + C x 2 , (5) 

substitution whereof in entails: 

Sxi + C x 2 = F . (6) 
Together with PJ, the latter yields: 

S - - t F(t/)x2(t/) dt> C - f F(t/)xi(t ' } dt> (7) 

J W[x 1 ,x 2 )(t') ' J W[x 1 ,x 2 )(t') [n 

where W^Xi, x 2 ]it) = Xi(t) x 2 (t) — x 2 (t) X\{t) . (8) 

This traditional way of introducing the method of variation of parameters (VOP) is peda- 
gogically flawed because it does not illustrate the full might and generality of this approach. 
What is important is that the initial equation, whose solution(s) is (are) assumed to be 
known, does not necessarily need to be linear. Moreover, the parameters to be varied 
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should not necessarily be the coefficients in the linear combination of solutions. Historically, 
Euler and Lagrange developed this approach in order to solve the nonlinear equation (127)1 . 
so the parameters to vary (the orbital Keplerian elements) were not coefficients of a linear 
combination of solutions to the homogeneous equation. Rather, these were the constants 
of the unperturbed motion (i.e., the quantities which were conserved in the homogeneous 
(2-body) case but no longer conserved in the inhomogeneous (N-body) case). 

A pretty evident aspect of this method is that any choice of the supplementary condition 
different from (j3J) will render the same solution y(x) , even though the expressions for the 
"constants" will differ from (jTJ). In other words, the solution y(x) remains invariant 

under a certain family of transformations of C\ t 2 • This internal symmetry is inherent in any 
problem that can be approached through the VOP method, including the N-body problem of 
celestial mechanics where the parameters Cj are the constants of the unperturbed motion. 
To illustrate this, let us consider a trivial example kindly offered to the author by William 
Newman. 

1.2 Newman's Example 

As a simplest instance, consider a forced harmonic oscillator 

H(x,p) = ^- + ^-+xF(t) 

that leads to the well known initial-condition problem 

x + x = F(t) , with x(0) and i(0) known. (9) 
As prescribed by the method of variation of parameters, we seek a solution in ansatz 

x = S(t) sin t + C(t) cos t . (10) 

which yields 

x = {S(t) sin t + C(t) cost} + S(t) cost- C(t) sin t . (11) 

The standard procedure implies that we put S(t) sint + C(t) cost = , in order to get rid 
of the ambiguity. The by-hand imposure of this equality is convenient but not necessarily 
required. Any other way of fixing the ambiguity, like for example, 

S(t) sin t + C(t) cost = $(t) (12) 

will be equally good. Then 

x = $ + S(t) cos t-C(t) sin t-S(t) sin t - C(t) cost , (13) 

whence 

x + x = $ + S(t) cos t-C{t) sin t (14) 

Thus one faces the system 

6 + S(t) cost -C(t) sin t = F(t) 

(15) 

S(t) sin t + C(t) cost = $(t) , 
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the first line being the equation of motion (obtained through combining (|14p with (jSj)), and 
the second line being identity (fT2l) Fl The system trivially resolves to 

S = F cost - 4 ($cost) 
dt v ; 

(16) 

C = -F sint + — ($sint) 
dt 

The function $ still remains arbitrary]! as can be easily seen either from the above 
derivation or from direct substitution of ([TBI) in (Q. Integration of (fTST) trivially yields 

S = J F cost' dt' - $cost + ci 

(17) 

C = - J F sin*' dt' + $sint + c 2 
inclusion whereof into ffTUl) entails 



x = S sin t + C cos t = — cost / F sin t' dt' + sint / F cost' dt' + Ci sint + c 2 cost (18) 

with the <3>-terms cancelled out. Analogy of this simple example with the osculating-element 
formalism is evident. Variable x considered in this example is analogous to the Cartesian 
coordinates and velocities in Section III, while parameters S and C are counterparts 
of the osculating elements. Equations ffTBT) are analogues to ( 16"61) - (J7JJ). These expressions 
for S and C through $ show that it is impossible (at least, in this particular example) 
to pick up such $ that the right-hand sides in both equations (Tl6|) simplify. Newman's 
example also thwarts one's hope to separate the time scales: as can be seen from (ITT)) , even 
when the perturbation F(t) is a slow function of time, the time evolution of the "osculating 
elements" S and C is determined by the "fast" time scale associated with sint and 
cost , i.e., with the solutions for the homogeneous equation. 

The above setting was much simpler than the N-body case where one is faced with 
nonlinearity and the parameters to be varied are not coefficients of a linear combination 
of solutions to a homogeneous equations. However, Newman's example well illustrates the 
origin of internal symmetry in the VOP context. This symmetry is a generic feature that 
shows itself in all implementations of the VOP approach, including the N-particle dynamics. 

Another good thing about Newman's example is that it illustrates the impossibility of 
time-scales' separation in the VOP method: in the right-hand side of (fTTj) the expressions 
standing under the integral contain both the free-oscillator and disturbing frequencies. No 
matter what gauge $ one chooses, at least one of these integrals will contain a mixture 
of frequencies. This way, if we have a fast oscillator subject to a slow force, we still cannot 
reduce its motion to a slow evolution of the amplitudes S(t) and C(t); instead, the time 
dependence of these amplitudes will contain fast terms (which is very counter-intuitive). 



^^This identity becomes a constraint if one fixes the functional form of function $ which thus far has 
remained arbitrary. 

2 Mind that this arbitrariness of <I> stays, no matter what the initial conditions are to be. Indeed, for 
fixed x(0) and i(0), the system C(0) = x(0) , $(0) + 5(0) = i(0) solves for S(0) and C(0) for an arbitrary 
choice of $(0). 
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2 Background 

Solar-System dynamics is, largely, variations of the old theme, the gravity law anticipated 
by Hook and derived from Kepler's laws by Newton: 

!fi 1 dUi 
E Grrij = rrii -— 1 , fy = p} - p { , i,j = 1,...,N , (19) 

rrii and p i being the masses and their inertial-frame-related positions, C7, being the overall 
potential acting on 

Ui = E Gm 3 — > (2°) 

the sign convention chosen as in the astronomical, not as in the physical literature. The 
equations of motion may be conveniently reformulated in terms of the relative locations 

fi = r is = ^ - p s , (21) 

p s standing for the position of Sun. The difference between 

% = £ G + g (22) 



rr. T ■ 



and 



& = E G^ + G^ (23) 

amounts to: 



h = V G^pl - E - G {m + ™ s) * = |5 (24) 



U being the new potential: 



U, = — + + R, , (25) 



with the disturbing function 



Ri = Y\ Gm, \\ - (26) 

m 3 14 p$ J 

singled out. 

Formulae ff24|) - ff25l) become trivial in the case of the two-body problem where only 
and m s are present. In this situation the disturbing function vanishes and the motion is, 
mathematically, equivalent to rotation about a nailed-down body of mass m ; + m s located 
at some fixed point O : 

r + 4~ = , f = p lanet - p s , fi = G (m + m s ) . (27) 
In here r = f\ = fj , because the subscript z runs through one value solely: i = 1 . 
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This setting permits exact analytical treatment that leads to the famous Newtonian 
result: the orbit is elliptic and has the gravitating centre in one of its foci. This enables 
a transition from Cartesian to Keplerian coordinates. For our further study this transition 
will be very important, so we shall recall it in detail. 

At any instant of time, the position r and velocity r of an orbiting body can be 
determined by its coordinates (x, y, z) and derivatives (x, y, z) in an inertial frame with 
origin located in point O where the mass irti + m s rests. The position of the orbital 
ellipse may be fully defined by the longitude of the node, Q ; the inclination, % ; and 
the argument of pericentre, to (instead of the latter, one can introduce the longitude of 
pericentre, Co = Q + u> ). The shape of the ellipse is parametrised by its eccentricity, 
e , and semimajor axis, a . Position of a point on the ellipse may be charachterised, for 
example, by the eccentric anomaly, E . As well known, 

E - e sin E = nt - B , (28) 

B being a constant of integration, and n being the mean motion defined as 

n = II 1 ' 2 a" 3 / 2 . (29) 

One can then introduce, following Kepler, the mean anomaly, M as 

M = E - e sin£ . (30) 

Let t D be the fiducial time. Then, by putting B = M Q + nt Q , we can introduce, instead 
of B , another integration constant, M Q . Hence, ([281 will read: 

M = M + n (t - t a ) , (31) 

the meaning of M Q being self-evident: it is the value of M at the reference epoch t a . So 
introduced the mean anomaly provides another parameterisation of the position of a planet 
on the ellipse. In the disturbed case the latter formula naturally becomes 



M 



M + [ n(t')dt' . (32) 



One more convenient parameter often employed in the literature is the mean longitude A 
defined by 

\ = u + M = n + uj + M. (33) 

Unless the inclination is zero, neither the longitude of the pericentre, uj , nor the mean 
longitude, A , is a true angle. They are sums of angles in two different planes that meet at 
the node. 

Planetary dynamics is based on application of the above, 2-body, formalism to the 
N-body case. Naively speaking, since the mutual disturbances of planets are very weak 
compared to the solar gravity, it seems natural to assume that the planets move along 
ellipses which are slowly evolving. Still, the weakness of perturbations is, by itself, a very 
shaky foundation for the varying-ellipse method. This so physically-evident circumstance 
has a good illustrative power but is of no help when the following questions arise: 

(1) To what degree of rigour can an orbit curve be modelled by a family of instantaneous 
ellipses having the Sun in one of their foci? Can this be performed exactly? 

(2) Is this representation of the curve by a family of ellipses unique? 
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These two questions will not seem anecdotal, if we recall that the concept of evolving 
instantaneous ellipses had been introduced into practice (and that major developments of the 
disturbing-function theory had been accomplished) long before Frenet and Serait developed 
the theory of curvesH (This historical paradox explains the reason why words "helicity" 
and "torsion" are still absent from the astronomy textbooks.) 

Fortunately, Lagrange, who authored the idea of instantaneous ellipses, fortified it with 
such powerful tools of calculus, that in this case they surpassed the theory of curves. More- 
over, these tools in no way relied on the weakness of the disturbances. Hence, Lagrange's 
treatment of the problem already contained an affirmative answer to the first question. 

Below we shall demonstrate that the answer to the second question is negative. Moreover, 
it turns out that the question calls into being a rich, though not new, mathematical structure. 
We shall show that the Lagrange system of equations for the instantaneous orbital elements 
possesses a hidden symmetry not visible with the naked eye. This symmetry is very similar 
to the gauge symmetry, one well known from electrodynamics. A careful analysis shows 
that the Lagrange system, as we know it, is written in some specific gauge: all trajectories 
constrained to some 9-dimensional submanifold in the 12-dimensional space constituted by 
the Keplerian elements and their time derivatives. 

Beside the possible practical relevance to orbit computation, the said symmetry unveils 
a fiber bundle structure hidden behind Lagrange's system of equations for the Keplerian 
elements. The symmetry is absent in the 2-body case, but comes into being in the N-body 
setting (N > 3) where each orbiting body follows a ellipse of varying shape, but the time 
evolution of the ellipse contains an inherent ambiguity. 

Here follows a crude illustration of this point. Imagine two coplanar ellipses sharing one 
focus. Let one ellipse slowly rotate within its plane, about the shared focus. Let the other 
ellipse rotate faster, also in its plane, in the same direction, and about that same shared 
focus. Suppose a planet is at one of the points of these ellipses' intersection. One observer 
may state that the planet is rapidly moving along a slowly rotating ellipse, while another 
observer may insist that the planet is slowly describing the fast-moving ellipse. Both de- 
scriptions will be equally legitimate, for there exists an infinite number of ways of dividing 
the actual motion of the planet into its motion along some orbit and simultaneous evolution 
of the orbit itself. Needless to say, the real, physical trajectory is unique. However, its 
description (parameterisation in terms of Kepler's elements) is not. A map between two dif- 
ferent (though physically-equivalent) sets of orbital elements is a symmetry transformation 
(a gauge transformation, in physicists' jargon). 

Lagrange never dwelled on that point. However, in his treatment he passingly introduced 
a convenient mathematical condition similar to (j3D, which removed the said ambiguity. This 
condition and possible alternatives to it will be the topic of the further sections. 

3 I am grateful to William Newman who drew my attention to this circumstance. 
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3 Keplerian coordinates in the 2-body and N-body 
problems: Osculating Elements vs Orbital Elements 

Although not widely recognised, the perturbation equations of celestial mechanics possess 
a gauge freedom. It is probable that this was already noticed by Euler and Lagrange in 
the middle of the XVIII century. However, although the existence of this freedom did not 
entirely escape attention^ its consequences have yet to be fully explored. 

Perhaps the easiest way to gain an appreciation of this freedom is to follow the derivation 
of the perturbation equations by application of the variation of parameters (VOP) technique 
as invented by Euler and Lagrange and shaped into its final form in Lagrange (1808, 1809, 
1810). Lagrange put the planetary equations in a closed form in which temporal derivatives 
of the orbital elements were expressed in terms of partial derivatives of the disturbing 
function with respect to the orbital elements. A closely related development is presented in 
the textbook by Brouwer and Clemence (1961). We shall start in the spirit of Lagrange but 
will soon deviate from it in two points. First, in this subsection we shall not assume that 
the disturbing force is conservative and that it depends upon the positions solely, but shall 
permit it to depend also upon velocities. Second, neither in this subsection nor further shall 
we impose the Lagrange constraint. 

Before addressing the N-particle case, Lagrange referred to the reduced 2-body problem, 

•• u r 

F + ^- = > 

(34) 

T = Tplanet T sun ; /J> = C(jTlp[ ane i TTlsun) 

whose generic solution, a Keplerian ellipse or a hyperbola, can be expressed, in some fixed 
Cartesian frame, as 



(35) 



or, shortly: 



X = 


fi (Ci, . 


■•) Cq, t) , 


X = 


9i (Ci, . 


.., Cq, t) 


y = 


h(C lt . 


.., Cq, t) , 


y = 




.., Cq, t) 


z = 


fa(C u . 


.., Cq, t) , 


z = 


9s (Ci,. 


.., Cq, t) 


— * 

r = 


f(Ci,.. 


■i Cq, t) , 


r = 


g(Ci,... 


, Cq, t) 



(36) 



the functions #j being, by definition, partial derivatives of /j with respect to the last argu- 
ment: 

C'=const 

Naturally, the general solution contains six adjustable constants, Ci, since the problem ([31 
is constituted by three, second order, differential equations. To find the explicit form of 



4 In 2002 a veteran of the Russian astronomy Yurii Batrakov mentioned to the author that back in 
his student years his lecturer at the St. Petersburg (then Leningrad) State University, Mihail Subbotin, 
emphasised in his course the possibility of making alternative choices. Interestingly, though, this fact was 
not developed in the published version of that lecture course (Subbotin 1968) where it is only passingly 
mentioned that the Lagrange constraint is imposed merely for convenience. 
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the dependence (E5D, one can employ an auxiliary set of Cartesian coordinates q , with an 
origin at the gravitating centre, and with the first two axes located in the plane of orbit. In 
terms of the true anomaly / , these coordinates will read: 

q 1 = r cos/ , q 2 = r sin/ , q 3 = . (38) 

In the two-body case, their time derivatives can be easily computed and expressed through 
the major semiaxis, a , the eccentricity, e , and the true anomaly, / : 

n a sin f n a (e + cos f ) , , 

* = • * = vi - e> ■ * = ■ (39) 

(true anomaly / itself being a function of a , e , and of the mean anomaly M = 
M Q + // o n dt , n = /i 1//2 a -3 / 2 ) . In the two-body setting, the inert ial- frame- related 
position and velocity will appear as: 

r = R(fl, i, cj) q(a, e, M D , t) , 

(40) 

r = R(0, i, a;) q(a, e, M Q , t) , 

R(fi, z, a;) being the matrix of rotation from the orbital-plane- related axes q to the 
fiducial frame (x, y, z) in which the vector r is defined. The rotation is parametrised by 
the three Euler angles: inclination, % ; the longitude of the node, Q ; and the argument of 
the pericentre, u . 

This is one possible form of the general solution (|36l) . It has been obtained under the 
convention that a particular ellipse is parametrised by the Lagrange set of orbital elements, 
Ci = (e, a, M a , to, Q, i) . A different functional form of the same solution is achieved in 
terms of the Delaunay set, Di = (L, G, H, uj, Q, M a ) . Still another possibility is to express 
the general solution through the initial conditions: then the constants (x Q , y Q , z Q , x Q , y Q , z a ) 
are the six parameters defining a particular orbit. The latter option is natural when the 
integration is carried out in Cartesian components, but is impractical otherwise. 

At this point it is irrelevant which particular set of the adjustable parameters is employed. 
Hence we shall leave, for a while, the solution in its most general form (1351 - [37|) and shall, 
following Lagrange (1808, 1809, 1810), employ it as an ansatz for a solution of the N-particle 
problem where the disturbing force acting at a particle is denoted by AF jf] 

? + 4 - = AF , (41) 

the "constants" now being time dependent: 

r = f (C^t),. ..,C 6 (t),t) , (42) 
and the functional form of f remaining the same as in fl36l) . Now the velocity 
dr <9f ^ <9f dCi _ ^ df dCi 



5 Our treatment covers disturbing forces AF that are arbitrary vector- valued functions of position and 
velocity. 
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will contain a new input, ^2(df / dCi)(dCi/ dt) , while the first term, g , will retain the same 
functional form as it used to have before. 

Substitution of f (C±(t), ...^C^it), t) into the perturbed equation of motion (Pill) gives 
birth to three independent scalar differential equations of the second order. These three 
equations contain one independent parameter, time, and six time-dependent variables Cj(i) 
whose evolution is to be determined. Evidently, this cannot be done in a single way because 
the number of variables exceeds, by three, the number of equations. This means that, 
though the physical orbit (comprised by the locus of points in the Cartesian frame and 
by the values of velocity in each of these points) is unique, its parameterisation in terms 
of the orbital elements is ambiguous. Lagrange, in his treatment, noticed that the system 
was underdefined, and decided to amend it with exactly three independent conditions which 
would make it solvable. 

Before moving on with the algebra, let us look into the mathematical nature of this 
ambiguity. A fixed Keplerian ellipse (!36|) . which is the solution to the 2-body problem (134"|) . 
gives birth to a time- dependent one-to-one (within one revolution period) mapping 

( d , ... , C 6 ) <— ( x(t) , y(t) , z(t) , x(t) , y(t) , z(t) ) . (44) 

In the N-body case, the new ansatz (|42p is incompatible with f[5o]) . This happens because 
now the time derivatives of coordinates G{ come into play in fj4*3j) . Hence, instead of (T4*4|) . 
one gets a time-dependent mapping between a 12-dimensional and a 6-dimensional spaces: 

( Ci(t) , ... , C 6 (t) , C l (t) , ... , C 6 (t) ) - ( x(t) , y(t) , z(t) , x{t) , y(t) , z(t) ) . (45) 

This brings up two new issues. One is the multiple scales: while the physical motion along 
an instantaneous ellipse parametrised by C{ is associated with the "fast time scale", the 
evolution of the osculating elements Ci(t) represents the "slow time scale". The quotation 
marks are necessary here, because in reality these time scales are inseparably connected, so 
that the "slow time scale" is not always slower than the "fast time scale" . What is important 
here is that, in general, ansatz (j42p gives birth to two separate time scales. The second 
important issue is that mapping (14 5 j) cannot be one-to-one. As already mentioned, the 
three equations of motion (1411) are insufficient for determination of six functions C\, ... Cg 
and, therefore, one has a freedom to impose, by hand, three extra constraints upon these 
functions and their derivatives^. 

6 In practice, the mean longitude A = A Q + J n(t) dt is often used instead of its fiducial-epoch value 
A . Similarly, those authors who prefer the mean anomaly to the mean longitude, often use M = M Q + 
J t n(t) dt instead of M a . While M a and A Q are orbital elements, the quantities M and A arc not. Still, 
the time-dependent change of variables from A D to A (or from M a to M ) is perfectly legitimate. Being 
manifestly time-dependent, this change of variables intertwines two different time scales: for example, M 
carries a "fast" time dependence through the upper limit of the integral in M = M Q + f. n(t) dt , and it 
also carries a "slow" time-dependence due to the adiabatic evolution of the osculating element M a . The 
same concerns A . 

7 A more accurate mathematical discussion of this freedom should be as follows. The dynamics, in the 
form of first-order differential equations for the orbital coordinates Ci{t) and their derivatives Hi(t) = 
Ci(t) , will include six evident first-order identities for these twelve functions: Hi(t) — dCi(t)/dt . Three 
more differential equations will be obtained by plugging r = f(Ci,...,Ce,t) into (j4Tj) . These equations 
will be of the second order in Cj(t) . However, in terms of both Cj(t) and Hi(t) these equations 
will be of the first order only. Altogether, we have nine first-order equations for twelve functions Ci(t) 
and Hi(t) . Hence, the problem is underdefined and permits three extra conditions to be imposed by 
hand. The arbitrariness of these conditions reveals the ambiguity of the representation of an orbit by 
instantaneous Keplerian ellipses. Mappings between different representations reveal an internal symmetry 
(and a symmetry group) underlying this formalism. 
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Though Lagrange did notice that the system was underdefined, he never elaborated on 
the underlying symmetry. He simply imposed three convenient extra conditions 



and went on, to derive (in this particular gauge, which is often called "Lagrange constraint") 
his celebrated system of equations for orbital elements. Now we can only speculate on why 
Lagrange did not bother to explore this ambiguity and its consequences. One possible 
explanation is that he did not have the concept of continuous groups and symmetries in his 
arsenal (though it is very probable that he knew the concept of discrete grourj^l). Another 
possibility is that Lagrange did not expect that exploration of this ambiguity would reveal 
any promising tools for astronomical calculations. 

Lagrange's choice of the supplementary constraints was motivated by both physical con- 
siderations and the desire to simplify calculations. Since, physically, the time-dependent 
set (C%(t), ...,C 6 (t)) can be interpreted as an instantaneous ellipse, in a bound-orbit case, 
or an instantaneous hyperbola, in a fly-by situation, Lagrange decided to make the instan- 
taneous orbital elements Cj osculating, i.e., he postulated that the instantaneous ellipse 
(or hyperbola) must always be tangential to the physical trajectory. This means that the 
physical trajectory defined by (Ci(t), C 6 (t)) must, at each instant of time, coincide with 
the unperturbed (two-body) orbit that the body would follow if perturbations were to cease 
instantaneously. This can be achieved only in case the velocities depend upon the elements, 
in the N-body problem, in the same manner as they did in the 2-body case. This, in turn, 
can be true only if one asserts that the extra condition (H6l) is fulfilled. That condition, 
also called Lagrange constraint, consists of three scalar equations which, together with the 
three equations of motion (HTi) . constitute a well-defined system of six equations with six 
variables (Ci (t) , . .. , Cq (t)) . 

For the reasons explained above, this constraint, though well motivated, remains, from 
the mathematical viewpoint, completely arbitrary: it considerably simplifies the calculations 
but does not influence the shape of the physical trajectory and the rate of motion along 
that curve. One could as well choose some different supplementary condition 

Z§/f = *<C .,,) , (47) 

$ now being an arbitrary function of time and parameters Ci E Substitution of (jlBjl by ( 1471) 
would leave the physical motion unchanged but would alter the subsequent mathematics 
and, most importantly, would eventually yield different solutions for the orbital elements. 
Such invariance of a physical theory under a change of parameterisation is an example of 
gauge symmetry (Efroimsky 2002). This freedom is a particular case of a more general type 



8 In his paper on solution, in radicals, of equations of degrees up to four, Reflexions sur la resolution 
algebrique des equations, dated by 1770, Lagrange performed permutations of roots. Even though he did 
not consider compositions of permutations, his technique reveals that, most likely, he was aware of, at least, 
the basic idea of discrete groups and symmetries. 

9 In principle, one may endow 3? also with dependence upon the parameters' time derivatives of all 
orders. This would yield higher-than-first-order time derivatives of the Ci in subsequent developments 
requiring additional initial conditions, beyond those on f and r, to be specified in order to close the system. 
We avoid this unnecessary complication by restricting <& to be a function of time and the Ci . 
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of invariance that emerges in applications of the VOP method to a vast variety of ODE 
problems. 

The importance of this gauge freedom is determined by two circumstances which parallel 
similar circumstances in electrodynamics. One consequence of the gauge invariance is its 
non-conservation in the course of orbit computation. This, purely numerical, phenomenon 
may be called gauge drift. Another relevant issue is that a clever choice of gauge often 
simplifies the solution of the equations of motion. In application to the theory of orbits, 
this means that a deliberate choice of non-osculating orbital elements (i.e., of a set C, 
obeying some condition (|47p different from (|46p ) can sometimes simplify the equations for 
these elements' evolution. As explained in Efroimsky & Goldreich (2003b), some earlier 
developments by Goldreich (1965), Brumberg et al (1971) can be interpreted as calculations 
in non-Lagrange gauges. 

The standard derivation of both Delaunay and Lagrange systems of planetary equations 
by the VOP method rests on the assumption that the Lagrange constraint (146!) is fulfilled. 
Both systems get altered under a different gauge choice. Derivation of the gauge-invariant 
(i.e., valid in an arbitrary gauge (1471) ) Lagrange system can be found in Efroimsky (2002). 
Derivation of the gauge-invariant Delaunay system is given in Efroimsky (2002). In Efroim- 
sky & Goldreich (2003b) both systems are amended with extra terms which emerge when 
the disturbance happens to depend not only upon the positions but also upon velocities. 
This situation is encountered when one is working in a non-inertial coordinate frame or 
when relativistic corrections to the Newtonian force are accounted for. 

Without going into excessive details, that may be found in (Efroimsky 2002), we give 
the essence of that derivation so that we can refer to it in one of the subsequent sections. 

^From (H3j) . it follows that the formula for the acceleration reads: 



d 2 f d£ _ dg_ dQ dj> dH ^ dj dCj d$> 

dt 2 dt ^ dQ dt dt d 2 t 4- dd dt dt ' 1 j 

Together with the equation of motion (T4T]) . it leads to: 

<9 2 f u f ^ dg dCi d<& ._. .rt. . 

^ + ^; + S4^ + ^ = ^ • r S |f| = |f| . (49) 

As f is, by definition, a Keplerian solution to the two-body problem (and, therefore, obeys 
the unperturbed equation ( l34l) ), the above formula becomes: 

This is the equation of disturbed motion, written in terms of the orbital elements. Together 
with constraint (1471) it constitutes a well-defined system of equations that can be solved 
with respect to dCi/dt for all z's. The easiest way of doing this is to employ the elegant 
technique offered by Lagrange: to multiply the equation of motion by df/dC n and to 
multiply the constraint by —dg/dC n . These operations yield the following equalities 

df ( d^dcA m_ f _ df_ d$ 



dC n \*f dCj dt j dC n dC n dt 
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and 

dC n dCj dt J dC n { ! 

summation whereof results in: 

V \ C C -\^ - —*F - — — - (53) 

\ [Cn CA dt ~ oc n Ab oc n dt dc n * ' (53) 

where the symbol [C n Cj] denotes the unperturbed (i.e., defined as in the two-body case) 
Lagrange brackets: 

[Cn Cj] = dc~ n dc 3 ~ dc^ dc~ n ' (54) 

Above we specified that $ is a function of time and the parameters Cj , but not of their 
derivatives. Under this restriction, (153]) may be conveniently rewritten as 

?H> + |;If)f "I * • w 

This expresses the most generic form of the gauge-invariant perturbation equations of ce- 
lestial mechanics, that follows from the VOP methodic 

We can obtain an immediate solution for the individual dCJdt in the Lagrange gauge by 
exploiting the well known expression for the inverse Lagrange-bracket, or Poisson-bracket, 
matrix that is offered in the literature for the two-body problem. The presence of the term 
proportional to d^/dCj on the lhs of (|65|) complicates the solution for the individual dCi/dt 
in an arbitrary gauge, but only to the extent of requiring the resolution of a set of six, 
simulataneous, linear, algebraic equations. 

All different choices of three (compatible and sufficient) gauge conditions expressed by 
the vector $ will lead to physically equivalent results. This equivalence means the fol- 
lowing. Suppose we solve the equations of motion for Ci,...^ , with some gauge condition 
$ enforced. This will give us the solution, C lr .. )6 (i) . If, though, we choose to integrate 

the equations of motion with another gauge $ enforced, then we shall arrive at a solution 
Ci(t) of a different functional form. Stated alternatively, in the first case the integration 
in the 12-dimensional space ( Ci,...^ , Hi,...,e ) will be restricted to one 9-dimensional 
time-dependent submanifold, whereas in the second case it will be restricted to another 
submanifold. Despite this, both solutions, Ci(t) and Ci(t) , will give, when substituted 
back in (H2"j) , the same orbit (x(t), y(t), z(t)) with the same velocities (x(t), y(t), z{t)) . 
This is a fiber-bundle-type structure, and it gives birth to a 1-to-l map of Cj(t) onto 
Ci(t) . This map is merely a reparameterisation. In physicists' parlance it will be called 
gauge transformation. All such reparameterisations constitute a group of symmetry, which 
would be called, by a physicist, gauge group. The real orbit is invariant under the reparam- 
eterisations which are permitted by the ambiguity of gauge-condition choice. This physical 
invariance implements itself, technically, as form-invariance of the expression (|42|) under the 
afore mentioned map. This is analogous to Maxwell's electrodynamics: the components 
x, y, and z of vector r, and their time derivatives, play the role of the physical fields E 



10 I am grateful to Peter Goldreich who recommended me to present the equations in this concise and so 
general form. 
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and B , while the Keplerian coordinates C\, C@ play the role of the 4-potential . This 
analogy can go even further 

4 The hidden symmetry of the Lagrange system 

If we impose, following Lagrange, the gauge condition (146]) . then the equation of motion 
( 165*1) will simplify: 

W.lf^AF . (56) 

In assumption of the disturbing force being dependent only upon positions and being con- 
servative!^ we may substitute in the above equation the disturbing force by the gradient of 
a (position-dependent) disturbing potential: 

which will result in: 

Expressions for the Lagrange brackets are known (Brouwer & Clemence 1961, p. 284), and 
their insertion into fl58l) equation will entail the well-known Lagrange system of planetary 
equation: 

da 2 dR 

(59) 
(60) 



(it na dM a 

de 1 -e 2 OR (1 - e 2 ) 1 ' 2 <9i? 
(it n a 2 e (9M n a 2 e 9a 



_ - cosz dR + (1-e 2 ) 1 / 2 

(it na 2 (l — e 2 ) 1 / 2 sin z di na?e de 



11 Suppose one is solving a problem of electromagnetic wave proliferation, in terms of the 4-potential 
in some fixed gauge. An analytic calculation will render the solution in that same gauge, while a 

numerical computation will furnish the solution in a slightly different gauge. This will happen because of 
numerical errors' accumulation. In other words, numerical integration will slightly deviate from the chosen 
submanifold. This effect may become especially noticeable in long-term orbit computations. Another 
relevant topic emerging in this context is comparison of two different solutions of the N-body problem: just 
as in the field theory, in order to compare solutions, it is necessary to make sure if they are written down 
in the same gauge. Otherwise, the difference between them may be, to some extent, not of a physical but 
merely of a gauge nature. 

12 I am grateful to Peter Goldreich who drew my attention to the fact that this assertion is right in inertial 
axes solely. Consider an observer who associates himself with a non-inertial system of reference (say, with 
a frame fixed on a precessing planet) and defines the orbital elements in this system. He will then have 
to take into account the non-inertial contribution to the disturbing function. This contribution, denoted 
in (Goldreich 1965) by Ri, will depend not only upon f but also upon r . In this situation, the terms 
dR/dC r in the Lagrange equations should be substituted by dR/dC r — (dr /dC r )(dR/df). The gauge 
approach to the problem of a satellite orbiting a wobbling planet will be addressed in one of our subsequent 
articles (Efroimsky & Goldreich 2003b) . In regard to the velocity-dependence of the disturbing function, we 
would also mention that a similar situation emerges in the relativistic mechanics, because the relativistic 
correction to the Newton law of gravity bears dependence upon the velocity. 
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di 
dJ 



na* 



cos i 
- e 2 ) 1 / 2 
dtt 
~dt 



dR 



dM n 



sin i olo n a 1 ( 1 — e 

1 dR 

na 2 (l — e 2 ) 1 / 2 sin i di 
1 - e 2 dR 2 dR 



J 1/2 



sin z 



(62) 
(63) 

(it n a 2 e de na da ^ 

If analytical integration of this system were possible, it would render a correct orbit, in 
the fixed gauge (1461) . A numerical integrator, however, may cause drift from the chosen 
submanifold f|46|) . Even if the drift is not steady, some deviation from the submanifold is 
unavoidable. One may wish <fr to be as close to zero as possible, but in reality $ will 
be some unknown function whose proximity to zero will be determined by the processor's 
error and by the number of integration steps. Even if we begin with (14*61) fulfilled exactly, 
the very first steps will give us such values of Cj that, being substituted into the lhs of 
( 1461) . they will give some new value of $ slightly different from zero. Thus, the Lagrange 
gauge will no longer be observed. 

In case we relax the Lagrange constraint and accept an arbitrary gauge (HTI) then, under 
the simplifying assertion ( 1571) . equation ( 1651) will become 



di d® 



dC n dCj 



dC± 
dt 



dR 

dCn. 



di d<& dg 

dC~ n ~dt ~ dC~ n 



(65) 



The La grange brackets depend exclusively on the functional form of x, y, z — /i,2,3 (Ci,...,6 > t) 
and #1,2,3 = dfi^/dt , and are independent from the gauge and from the time evolution of 
Cj. Hence the gauge- invariant generalisation of the Lagrange system will emerge: 



da 
~d~b 



2 

na 



dR 
~dM 



di d<& 



dM D dM Q dt 



(66) 



de 
dt 



na* e 



dR 
dMl 



_dg_ 

dM n 



di d<f> 



dM dt 



3 2U/2 



n a 2 e 



dR 
da 



da 



di d<f> 

da dt 



(67) 



du 
~dl 



— COS I 



na* 



2U/2 



e 2 ) 



sin i 



dR ^ dg di d<f> 

di di di dt 



3 2U/2 



na 2 e 



+ 



dR ^ dg di d& 

de de de dt 



(6* 



di 
dt 



cos % 



n a 2 (1 — e 2 ) 1 / 2 sin i 



dR ^ dg di d& 

dcu duj dui dt 



dR - dg 
dti ~ dtt 



na 2 (1 — e 2 ) 1 / 2 sin i 



di_ d$ 

dVt ~dt 



(69) 
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dn 1 

dt na 2 (l — e 2 ) 1 / 2 sin i 

dM Q __ 1 - e 2 
dt na 2 e 

These gauge-invariant equations reveal the potential possibility of simplification of orbit 
integration. One can deliberately choose gauges different from (j4~6]) . In principle, it is 
possible to pick up the gauge so as to nullify the right-hand sides in three out of six equations 
(16"6l - I7T|) . This possibility is worth probing (we know from electrodynamics that a clever 
choice of gauge considerably simplifies solution of the equations of motion). Slabinski (2003) 
recently offered an example of gauge that nullifies the rhs of equation (1691) . 

— * 

Another tempting possibility may be to pick up the gauge so that the $-terms in (1661 
- ITTj) fully compensate the short-period terms of the disturbing functions, leaving only the 
secular and resonant ones. However, the afore discussed Newman's example strongly speaks 
against this possibility, at least in the case of bound orbits. The case of flybys seems to be 
more favourable, because in that case we do not have two different time scales. Hence in that 

— * 

case a choice of some nonvanishing <fr may, potentially, lead to an even larger simplification 
of calculations. We shall address this matter in a separate work. 

— *■ 

Finally, we would point out that, in case the gauge <1> depends not only upon time but 
also upon the "constants" Cj , the right-hand sides of (1661) - (1711) implicitly contain time 
derivatives of . Then, in order to continue with analytic developments, it is necessary to 
transfer such terms to the left-hand side, as we did in (1651) . This alteration will be dwelled 
upon in (Efroimsky & Goldreich 2003a,b). 



dR 

di 



- $ 



dg <9f r/<& 



di 



di dt 



(70) 



dR 

de 



de 



di d<S> 

de dt 



2 

na 



dR 
da 



dg 

da 



df d<S> 

da dt 



(71) 



5 Delaunay's variables 

As well known (Brouwer & Clemence 1961, p. 290), it is possible to choose as the parameters 
Ci not the six Keplerian elements Ci = ( e , a , M D , u , Q , i) but the set Di = 
( L , G , H , M Q , u , Q ) , new variables L , G , and H being defined as 

L^/i^a 1 / 2 , G = //V/ 2 (l - e 2 ) 1/2 , H = /i 1/2 a 1 / 2 (l - e 2 ) V ' cos i ,(72) 

where fi = G{m sun + m p i anet ) . 

The advantage of these, Delaunay, variables lies in the diagonality of their Lagrange- 
bracket matrix. Inversion thereof yields, similarly to (I66p - ( 171 p . the so-called Delaunay 
system: 



dL 


dR 


dM 


dR 


~dt ~ 


dM ' 


dt 


~~ dL 


dG 


dR 


du 


dR 


~dt 


dcu 


~dt ~ 


~ dG 


rill 


dR 


dn 


dR 


~dt 




~dt ~ 


~ ~dH 



(73) 



provided these parameters obey the Lagrange- type gauge condition analogous to (jlBl) : 

E^^ = 0- (74) 
^ dD f dt 1 ' 
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where, similarly to f[3"o1 - I3T1) . 



~ op 

f = f (D lr .., 6 , t) , 'r = i(D ly .., 6 ,t) , % = i , (75) 

and the tilde symbol emphasises that the functional dependencies of the position and velocity 
upon Di differ from their dependencies upon Cj . We, thus, must keep in mind that the 
system of equations for the Delaunay elements only pretends to exist in a 6-dimensional 
phase space. In reality, it lives on a 9-dimensional submanifold (1711) of a 12-dimensional 
manifold spanned by the Delaunay elements and their time derivatives. In the case of 
analytical calculations this, of course, makes no difference. Not in the case of numerical 
computation, though. 

If instead of the gauge condition (IT4I) we impose some alternative gauge 

^ df dD,, . . 

i L 

the generalised Delaunay-type equations will read@ 

dL dR | di df dM Q dR |<9| df d® 

~dt ~ ~dM ~ d~M Q ~ dM~ ~dt ' ~dT ~ dL + dL + dL ~dt ' 

dG dR ^di df d$ du dR ^dg df d® 

= — — = + — + (77) 

dt dcu du dcu dt ' dt dG dG dG dt ' 1 ; 

dH dR | di df d® dtt dR | di df d$ 

~dT ~ dQ ~ dQ ~ M ~dt ' ~dT ~ ~ dH + dH + dH ~dT 1 

and the $ terms should not be ignored, because they account for the trajectory's deviation 
from the submanifold f!74|) of the ambient 12-dimensional space (Di,...$ , Di,...,oj ■ 

The meaning of f and g in the above formulae is different than that of f and g 
in Sections 3 and 4. In those sections f and g denoted the functional dependencies 
(l3"6l) of x , y , z and x , y , z upon parameters = ( e , a , M G , a; , ft , i ) . Here, 

f and g stand for the dependencies of x , y , z and x , y , z upon the different set 
Di = (L , G , H , M , u; , f2 ) . Despite the different functional forms, the values of f 
and g coincide with those of f and g : 

?Pi,.., 6 ) = f f = f(Ci,... j6 ) and i(A,...,6) = f = g(d 6 ) . (78) 

Similarly, $ (-Di,...,6) an d $ (Ci,...^) are different functional dependencies. It is, though, 
easy to show (using the differentiation chain rule) that their values do coincide: 

i(D lt ^ 6 ) = $(^,...,6) (79) 



13 Similarly to our comment in the end of Section 4 we would mention that if the gauge $ depends not 
only on time but also on the parameters Cj then the right-hand sides of (|77p contain time derivatives of 
Q . A further analytic treatment of (|77p will then demand that we transfer those terms to the left-hand 
side, as in ([55)1 . This will be discussed in (Efroimsky & Goldreich 2003a, b). 
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which is analogous to the covariance of Lorentz gauge in electrodynamics. This means that, 
for example, analytical calculations carried out by means of the Lagrange system (1591 - [6^1) 
are indeed equivalent to those performed by means of the Delaunay system (1731) . because 

imposure of the Lagrange gauge $ = is equivalent to imposure of $ = . 

Can one make a similar statement about numerical integrations? This question is non- 
trivial. In order to tackle it, we should recall that in the computer calculations the Lagrange 
condition $ = cannot be imposed exactly, for the numerical error will generate some 
nonzero $ . In other words, the orbit will never be perfectly constrained to the submanifold 
$ = . Thereby, some nonzero $ will, effectively, emerge in (I661- I7TT) . Similarly, a small 

nonzero will, effectively, appear in (TM1) . and the Delaunay system will no longer be 
canonical. In other words, we get not just an error in integration of the canonical system, 
but we get an error that drives the system of equation away from canonicity. This effect 
is not new: it is well known that not every numerical method preserves the Hamiltonian 

structure. Therefore, the unavoidable emergence of a nonzero numerical-error-caused $ in 
the Delaunay system may, potentially, be a hazard. This topic needs further investigation. 

The gauge-invariant Delaunay-type system (I77|) is no longer Hamiltonian, unless the 
gauge is that of Lagrange, $ = . It is possible to demonstrate that by introducing a 
velocity-dependent disturbing force (which can be done, for example, by defining the orbital 
elements in a non-inertial reference frame) we further alter the Lagrange-type and Delaunay- 
type gauge-invariant systems. The gauge-invariant Delaunay-type system will remain non- 
canonical, except in one special gauge which is called in (Efroimsky & Goldreich 2003a) 
"the generalised Lagrange gauge" and which simplifies to the regular Lagrange gauge when 
the velocity dependence is removed. Applications of this construction will be explained in 
Efroimsky & Goldreich (2003b). 

To draw to a close, we would mention that in most books the Lagrange and Delaunay 
systems of equations are derived not by a straightforward application of the VOP method 
but through the medium of Hamilton- Jacobi technique (Plummer 1918, Smart 1953, Pollard 
1966, Kovalevsky 1967, Stiefel and Scheifele 1971). At first glance it may seem that such a 
derivation does not exploit any sort of supplementary constraint and is, therefore, spared 
of the gauge invariance. In our upcoming publication Efroimsky & Goldreich (2003a) we 
shall explain that the Lagrange constraint is implicitly imposed within the Hamilton- Jacobi 
treatment of this problem. Another method of derivation, different from both the VOP and 
the Hamilton- Jacobi methods, was offered in the book by Kaula (1968) and, once again, it 
may seem that his development is devoid of whatever extra conditions. In the next section 
we explain how Kaula's method, too, tacitly exploits the Lagrange constraint. 



6 Brouwer and Kaula 

As we saw above, the gauge-invariant generalisation of the Delaunay system is no longer 
canonical. In celestial mechanics, equations of type 

r = ~dp + 

(80) 

P = ~ ~dr + ^' P ' 
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were introduced by Dirk Brouwer who used them to include the atmospheric-drag forces 
into the canonical picture. The quantities X and P are called "canonical forces". The 
appropriate formalism is comprehensively explained in Stiefel & Scheifele (1971). We see 
now that such "forces" can emerge not only due to real physical interactions but also as 
purely mathematical artifacts called into being by a nontrivial gauge choice. This in no 
sense means that such terms may be omitted. Their relevance to the stability of numerical 
methods will be discussed in Murison & Efroimsky (2003). 

Above we saw that for a disturbance of type (j57]l there exists one special gauge (the 
Lagrange one) that makes the equations for Delaunay's variables canonical. Still, in the 
case of an arbitrary gauge the symplectic structure is destroyed!^! To better understand 
the anatomy of this phenomenon, let us explore an alternative derivation of Delaunay's 
equations, one based on a direct change of variables. This method is akin to the VOP 
technique, but has a better illustrative power. To simplify things, we shall assume that the 
disturbing force depends on the positions solely. Following Kaula (1968), we shall begin 
with the canonical equations for the position and velocity in some inertial Cartesian frame: 

dr ^ 

dt = r ' 

(81) 

d'r _ dU 
dt dr 

A solution to these will depend upon seven quantities So = t, Si, ... , Se ■ One of these will 
be time, another six will be the parameters which can be devoid of or imparted with a time 
dependence of their own. This paves our way to the new, modified, system of equations: 



dxi dSh 



dS k dt 
dfi dS k dU 



(82) 



dSk dt dri 

summation over repeating indices, from through 6 , being implied. After we multiply 
the upper and the lower equations by — dri/dSi and dri/dSi , correspondingly, we should 
sum up both lines and arrive to: 

[Si > Sk] -dr = ds l { u - 2™) ' [Si > Sk] " ait di - ait di ■ (83) 

What remains now is to calculate the Lagrange-bracket matrix [Si , S&] and to invert the 
result. This will yield the desired expressions for dS^/dt . The expression for k = 
will be merely an identity 1 = 1, because So = t . The other six will be the planetary 
equations. In case the parameters Si, ... , Sg make a Delaunay set, these equations will be 
expected to have a Hamiltonian form. Have they really? To check this, one has to calculate 
the appropriate Lagrange brackets. Or to accurately examine how these are calculated in the 
literature. The most direct way of computing the brackets is to differentiate formulae (1401) 



14 It will be proven in our next paper (Efroimsky & Goldreich 2003a) that in the more general case 
of a disturbing force dependent on positions and velocities the situation is similar: the equations for the 
Delaunay elements will become symplectic in a so-called generalised Lagrange gauge. 
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with respect to the orbital elements. This will give the well-known time-independent expres- 
sions for the Lagrange brackets. These often-quoted standard expressions will, though, be 
valid in the two-body case solely. The latter circumstance is missing in (Kaula 1968) where 
the author implies that all will work equally well in the case of N bodies. In reality, two 
types of additional inputs will show themselves in the above formula for r in the N-body 
problem: First of all, extra items will appear in the expression for q: 

l, <9q dqda dqde dq dM Q 

q ~ ~dt + dadt + de~~dl + dM~ ~dT ' ^ ^ 

Beside this, the contribution Rq will now enter the expression for x . Here 

• dUdil dUduj dUdi 

mUt + IkJ ~dt + ~didt ' ( ' 

Altogether, these will result in the following expression for the velocity: 

'r = R(fi, i, u) c[(a, e, M) + Rq 

R dq <9(Rq) da d(Rq) de d(Rq) dM Q fl(Rq) dtl g(Rq) du g(Rq) di 

dt da dt de dt dM Q dt dVt dt doj dt di dt 

- R§ + * . (86) 

In Kaula (1968) the latter term is ignored, so that the author implicitly imposes the Lagrange 
gauge. In this gauge, differentiation of the above expression with respect to the elements 
will indeed entail the standard time-independent expressions for the Lagrange brackets. In 
a non-Lagrange gauge, though, the situation will be different, and we shall not end up with 
a canonical Delaunay system. Instead, we shall arrive to the non-canonical gauge-invariant 
Delaunay-type system (!77|) . 

Very similarly, the Lagrange constraint enters all the methods by which the Gauss system 
of equations is derived in the literature, and the imposure of this constraint is camouflaged 
in a manner similar to what we saw above in regard to the Delaunay equations. We shall not 
engage in a comprehensive discussion of this issue, but shall rather provide a typical example. 
In section 11.5 of (Danby 1988) an auxiliary vector h is introduced as a unit vector aimed 
along the instantaneous orbital momentum of the body, relative to the gravitating centre. 
Then h gets resolved along the inertial axes (x, ,y, z) , where x points toward the 
vernal equinox and z toward the north of the ecliptic pole: 

h = x sin Q sin i — y cos Q sin i + z cos % . 

This expression is certainly correct in the two-body case. It remains valid also in the N-body 
problem, only if the orbital elements Q and i are osculating, i.e., only if the instantaneous 
inertial velocity is tangential to the ellipse parametrised by the Keplerian set that includes 
these fl and i . Thus, the Lagrange constraint is implied. 
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7 Is it worth it? 



At this point one may ask if it is at all worth taking the nonzero $ into account in the 
Delaunay equations. After all, one can simply restrict himself to the 6- dimensional phase 
space defined by Di , and postulate that the six unwanted extra dimensions Di do not 

exist (i.e., postulate that $ = ). This, of course, can be done, but only at some cost: 
a certain type of accumulating numerical error will be ignored (not eliminated), and it will 
keep accumulating. As explained in the end of the previous section, the overall integration 
error of a Hamiltonian system consists of an error that leaves the system canonical (like, for 
example, an error in calculation of R in (1771) ) and an error that drives the system away 
from its canonicity (like the error reflected in the accumulated nonzero value of $) . The <fr 
terms in (I77j) play the role of correctors: they fully compensate for the errors of the second 
type (i.e., for what in numerical electrodynamics is called gauge shift). 

Similarly, in the case of Lagrange system, one may enquire if it is worth introducing the 
12-dimensional space spanned by the orbital elements Cj and their time derivatives Hi = Ci. 
Why not simply consider a trajectory in the 6-dimensional space of C$ and assume that 
the Lagrange gauge is fixed exactly? Indeed, if we are solving the problem r = f (r) , is it 
worth introducing an extra entity v = r and considering the orbit integration in the space 
of a larger dimension, spanned by the components of both r and v ? Will this new entity 
v add anything? The answer to this question will be affirmative if we take into account 
the fact that a trajectory is not merely a locus of points visited by the body: the notion 
of trajectory also comprises the rate at which the body was travelling. Appropriately, the 
accumulated numerical error will consist of two parts: distortions of the orbit shape, and 
distortions in the time-dependence of the speed at which the orbit was described. This 
explains why the events are taking place not just in the space of orbital elements but in the 
larger space of the elements and their time derivatives. This issue is best of all explained 
by Hagihara (1970). In subsection 1.6 of the first volume of his treatise he contrasts the 
cases of exact and apparent equivalences of dynamical systems. (The latter case is that of 
the orbital curves being geometrically, not dynamically, identical.) 

Still, there is more to it, because a convenient choice of gauge may simplify the solution 
of the equations of motion (Slabinski 2003; Efroimsky & Goldreich 2003b). 

8 Conclusions 

We have demonstrated a previously unrecognised aspect of the Lagrange and Delaunay sys- 
tems of planetary equations. Due to the Lagrange gauge condition (116]) . the motion is, in 
both cases, constrained to a 9-dimensional submanifold of the ambient 12-dimensional space 
spanned by the orbital elements and their time derivatives. Similarly to the field theory, 
the choice of gauge is vastly ambiguous and reveals a hidden symmetry (and an appropriate 
symmetry group) inherent in the description of the N-body problem in terms of the orbital 
elements. Just as a choice of a particular gauge simplifies solution of the equations of motion 
in electrodynamics, an alternative (to that of Lagrange) choice of gauge can simplify orbit 
calculations. We have written down the generalised Lagrange-type ([66] -HQ and Delaunay- 
type (1771) equations in their gauge- invariant, form (with no specific gauge imposed). We have 
pointed out that neither the Lagrange gauge (jl6"]) nor any other constraint is exactly pre- 
served in the course of numerical computation. This may be a source of numerical instability. 
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